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The dynamic structure factor, vorticity and entropy density dynamic correlation functions are 
measured for Stochastic Rotation Dynamics (SRD), a particle based algorithm for fluctuating flu- 
ids. This allows us to obtain unbiased values for the longitudinal transport coefBcients such as 
thermal diffusivity and bulk viscosity. The results are in good agreement with earlier numerical and 
theoretical results, and it is shown for the flrst time that the bulk viscosity is indeed zero for this 
algorithm. In addition, corrections to the self-diffusion coefficient and shear viscosity arising from 
the breakdown of the molecular chaos approximation at small mean free paths are analyzed. In 
addition to deriving the form of the leading correlation corrections to these transport coefficients, 
the probabilities that two and three particles remain collision partners for consecutive time steps 
are derived analytically in the limit of small mean free path. The results of this paper verify that we 
have an excellent understanding of the SRD algorithm at the kinetic level and that analytic expres- 
sions for the transport coefficients derived elsewhere do indeed provide a very accurate description 
of the SRD fluid. 

PACS number(s): 47.11. -fj, 05.40.-a, 02.70.Ns 

I. INTRODUCTION 

Several years ago, Malevanets and Kapral (ij, \^ derived a simple and appealing algorithm — often called Stochastic 
Rotation Dynamics (SRD) or Multi-Particle Collision Dynamics (MPCD) — for the mesoscale modeling of fluctuating 
fluids. SRD is a particle based simulation technique with simple discrete time dynamics consisting of consecutive 
streaming and collision steps. It shares many features with Bird's Direct Simulation Monte Carlo (DSMC) algo- 
rithm 0, but uses more efficient multi-particle collisions to exchange momentum between the particles. Since there 
is a Boltzmann _ff-theorem for the SRD algorithm, and the particle number, momentum, and energy are locally 
conserved, the correct hydrodynamic behavior is guaranteed at large length and time scales. The algorithm therefore 
provides a convenient computational tool for solving the underlying thermo-hydrodynamic equations by providing 
a "hydrodynamic heat bath" which incorporates thermal fluctuations and provides the correct hydrodynamic inter- 
actions between embedded particles or polymers. An important advantage of SRD is that its simplified dynamics 
has enabled the analytical calculation of the transport coefficients and made it possible to obtain a rather complete 
theoretical understanding of the time-dependent correlation functions, the relaxation to equilibrium 0, , and the 
behavior in shear flow, including shear thinning at high shear rates [g. Because the algorithm correctly includes 
long-ranged hydrodynamic interactions and Brownian fluctuations — both of which are generally required for a proper 
statistical treatment of th e dy namics of mesoscopic suspended particles — it has been used to study the behavior of 
polymers 0, S ; colloids [lOl ITll | , including sedimentation 0, 0, 0| , and vesicles in flow 0, ^ . 

In its original form I', the SRD algorithm was not Galilean invariant at low temperatures, where the mean 
free path. A, is smaller than the cell size a. However, it was shown 0, 0| that Galilean invariance can be restored 
by introducing a random shift of the computational grid before every collision. In addition to restoring Galilean 
invariance, this grid shifting procedure accelerates momentum transfer between cells and leads to a coUisional contri- 
bution to the transport coefhcients. Two approaches have been used to analyze the resulting algorithm and determine 
the shear viscosity and thermal diffusivity. In Refs. and 0, a non-equilibrium kinetic approach is used to de- 
rive the transport coefficients. In Refs. ^ discrete-time projection operator technique was utilized to obtain 
Green-Kubo relations for the model's tr ansp ort coefficients, and explicit expressions for the transport coefficients were 
derived in accompanying papers '2l|'22'|. The two approaches are complementary and, for the most part, agree 
in their conclusions. The first is rather straightforward and intuitively appealing, but makes several assumptions 
which are not easily verified. The projection operator approach justifies in detail several assumptions used in the 
non-equilibrium calculations of Refs. 5] and 18]; it can also be used to analyze the transport coefficients of the 
longitudinal modes, namely the bulk viscosity and thermal diffusivity, which are hard to calculate analytically in the 
non-equilibrium approach [^. Note, in particular, that the collisional contribution to the thermal conductivity has 
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not yet been determined using the non-equilibrium methods. 

In spite of some claims to the contrary (5||, both approaches yield the same results for the transport coefficients. 
Table ^ contains a summary of both the collisional and kinetic contributions to the transport coefficients, including 
references to the original source of the results; the table caption also provides a brief synapsis of some misprints 
in published articles which may lead to confusion. Simulation results have generally been in good agreement with 
these predictions. This isparticularly true for the shear viscosit^where the results of equilibrium measurements of 
vorticity fluctuations [3jl7[ and the Green-Kubo relations |l9l |20| are in excellent agreement with non-equilibrium 
measurements in shear [5|,ll8] and Poiseuille 23] flow. The situation with the transport coefficients of the longitudinal 
modes, namely the bulk viscosity and thermal diffusivity, is somewhat less clear. The only reliable measurements of the 
thermal diffusivity have entailed equilibrium measurements of the corresponding Green-Kubo relation I20I I21I |23 | 
and a non-equilibrium measurement obtained by setting up a temperature gradient and measuring the resulting 
energy flux in a regime where the collisional contribution is negligible ||5j]. To our knowledge, there has been no direct 
verification of the prediction that the bulk viscosity is zero for SRD. 

In this paper we take an alternative approach based on an analysis of the equilibrium fluctuations of the hydro- 
dynamic modes to directly measure the shear and bulk viscosities and thermal diffusivity. Starting with an analysis 
of vorticity fluctuations to determine the shear viscosity, measurements of the dynamic structure factor are then 
used to deduce the values of the speed of sound, the thermal diffusivity, and the bulk viscosity. Measurements of 
the temporal behavior of the entropy correlations are also used to obtain a direct independent measurement of the 
thermal diffusivity. To our knowledge, this is the first quantitative measurement of the dynamic structure factor for 
SRD. An earlier measurement by Inoue et al. p3 | lead to unphysical results in the large frequency limit, and could 
not be used to determine the transport coefficients. 

The results of these measurements verify directly that the bulk viscosity is indeed zero for this algorithm. In addition, 
in agreement with earlier work UtI ITsl HqI f20l l2l| , results for the shear viscosity and the thermal diffusivity are 
in excellent agreement with the theoretical predictions presented in Tabled for a wide range of particle densities and 
mean free paths. However, as noted originally in Ref. jl^, and discussed in more detail in Ref. correlations 
between particles occupying the same collision cell at different time steps lead to an enhanced kinetic contribution 
to the transport coefficients. This breakdown of the molecular chaos approximation becomes pronounced at small 
mean free path. A, since particles do not travel far between collisions and tend to repeatedly have the same collision 
partners. For most transport coefficients, this additional contribution to the transport coefficients is masked by the 
collisional contribution, which dominates in the small mean free path regime. The affect is particularly pronounced, 
however, for the self-diffusion coefficient, for which there is no collisional contribution. Indeed, RipoU et al. p5l l2^ 
have observed that the self diffusion coefficient is significantly larger than the theoretical prediction of Ref. [l7ll2Cll | for 
small A/a. RipoU et al. provided a semi-analytical description of this behavior in which they determined numerically 
the number of particles sharing the same cell as a function of time. In this paper, we provide a detailed discussion 
of the leading correlation corrections to the kinetic contribution of both the shear viscosity and the self-diffusion 
coefficient and determine analytically the probability that two and three particles are in the same collision cells for 
consecutive time steps. While our results for the self-diffusion coefficient are in general agreement with those of Ref. 
|26j |. there seem to be several misprints in ^J, making a detailed comparison difficult. 

The remainder of the paper is organized as follows. After a brief summary of the SRD algorithm in Sec. II, 
the hydrodynamic equations of a simple liquid are reviewed and the correct form of the constitutive equations are 
discussed in Sec. III. The consequences of the fact that angular momentum is not conserved in the SRD algorithm 
are summarized, and the correct form of the viscous stress tensor is discussed. In particular, it is emphasized that in 
two dimensions, there is no difference between the viscous stress tensor of a simple isotropic fluid and an SRD fluid. 
The slight differences in three dimensions leave the form of the Navier-Stokes unchanged, with only a reinterpretation 
of the coefficient of sound attenuation. Sec. IV contains a fairly detailed derivation of the dynamic correlation 
functions in a simple liquid. The discussion follows rather closely that of Ref. jl^, but is included because several 
aspects of the derivation are a bit subtle and are generally not addressed in the literature. Explicit expressions for 
the vorticity, density, and entropy density dynamic correlations functions are presented. In Sec. V, these results 
are used to determine the shear and bulk viscosities and the thermal diffusivity. The agreement with the theoretical 
predictions summarized in Table is excellent. Sec. VI contains a detailed discussion of the consequence of the 
breakdown of the molecular chaos approximation at short mean free paths. While correlation effects do not change 
the collisional contributions to the transport coefficients .21.] . they do dramatically increase the amplitude of the 
kinetic contribution. In addition to deriving the form of the leading correlation contributions to the shear viscosity 
and the self-diffusion coefficient, the probabilities that two (^2) and three (ps) particles are in the same collision 
cells at consecutive time steps are derived analytically in the limit A/a — > 0. More generally, simulation results are 
used to show that p2 is solely a function of A/a. In the case of the shear viscosity, it is shown that inclusion of the 
leading correlation corrections yields results in surprisingly good agreement with measurements of the viscous stress 
correlations. For the self-diffusion coefficient, however, the correlation corrections for larger time intervals are large. 
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and dramatically increase the measured value of the self-diffusion coefficient for mean free paths smaller than the 
cell size. It is important to note that the correlation corrections considered here — arising from particles which collide 
with the same particles in consecutive time steps — are similar to those which occur in dense fluids interacting through 
soft potentials, and should therefore be interpreted as a "potential" or "coUisional" contribution to the velocity or 
stress correlation functions rather than a precursor of the power-law (long-time) tails observable at longer times. We 
believe that it is important to distinguish between these two effects, since long time tails are also visible at large mean 
free paths where these corrections are negligibly small. Although the same approach can be used to calculate these 
contributions, the corresponding probabilities are much harder to estimate. 

The results of this paper verify that we have an excellent understanding of the SRD algorithm at the kinetic level 
and that — with the exception of the self-diffusion coefficient — the analytic expressions for the transport coefficients 
given in Table ^ do indeed provide a very accurate description of the SRD fluid. Furthermore, the analysis of the 
dynamical structure factor and the dynamic entropy density correlation function verify directly that the algorithm 
satisfies the fluctuation-dissipation theorem. While this is to be expected for the current algorithm, which satisfies the 
necessary semi-detailed balance conditions 0, ^3 > verification studies of this type will be important for generalizations 
of the algorithm which model excluded volume effects through the use of biased multi-particle collision rules which 
depend on the local velocities and densities j28i i22|| . 

II. MODEL 

In the SRD algorithm, the fluid is modeled by particles with continuous spatial coordinates ri{t) and velocities 
Vi(t). The system is coarse-grained into the cells of a regular lattice with no restriction on the number of particles in a 
cell. The evolution of the system consists of two steps: streaming and collision. In the streaming step, the coordinate 
of each particle is incremented by its displacement during the time step, r. Collisions are modeled by a simultaneous 
stochastic rotation of the relative velocities of every particle in each cell. As discussed in Refs. Q and |3, a random 
shift of the particle coordinates (or, equivalently, the cell grid) before the collision step is required to ensure Galilean 
invariance. All particles are shifted by the same random vector with components in the interval [—a/2, a/2] before 
the collision step. Particles are then shifted back to their original positions after the collision. If we denote the cell 
coordinate of the shifted particle i by the algorithm is summarized in the equations 

r,{t + T) = r,{t)+Tv,{t) (1) 
v,(i + r) = nm + r)] + + r)] • {v,(t) - u[|^(i + r)]}, (2) 

where uj(^^) denotes a stochastic rotation matrix, and u(^|) = -p- is the mean velocity of the particles in 

cell All particles in the cell are subject to the same rotation, but the rotations in different cells are statistically 
independent. There is a great deal of freedom in how the rotation step is implemented, and any stochastic rotation 
matrix consistent with detailed balance can be used. In two dimensions, the stochastic rotation matrix, a;, is typically 
taken to be a rotation by an angle ±a, with probability 1/2 (see Refs. 0, ^3 01)- In three dimensions, one can 
perform rotations by an angle a about a randomly chosen direction, where all orientations of the random axis occur 
with equal probability (Model A in Ref. ^2^). 

III. HYDRODYNAMICS AND TRANSPORT COEFFICIENTS 

There is a hydrodynamic mode associated with each conserved density in a fluid. For a simple liquid, the conserved 
quantities are the particle mass density, p{r,t), the momentum density, g;(r,t), and the energy density, e(r, i), and 
the corresponding microscopic conservation laws are 

dtpiv,t)+dpgf3ir,t)^0, (3) 

9tg„(r,i) + 5/3f„0(r,t) =0, (4) 

and 

dte{r,t) + df3X0{v,t) = 0, (5) 

where TQ^(r, t) are the Cartesian components of the microscopic stress tensor and Xa(r, t) is the a-component of the 
microscopic energy current density. While Eqs. Q-© are microscopically exact, macroscopic constitutive relations 
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are required to close the system of equations. The constitutive equations relate local non-equilibrium averages of g(r, t) , 
'ra^(r,i), and Xa{r,t) to the local hydrodynamic variables, p{r,t) = (/5(r, i)), g(r,t) = (g ( r, t ) ), e{r,t) = (e(r, i)), and 
their gradients. For a simple isotropic liquid, the constitutive relations have the form [23, 123 



(ga(r,t)} = 9a{r,t), (6) 

{Tai3{r,t)) = p{r,t)Saf) - (Jai3{r,t) (7) 

= p{r,t)Sap - '^[dagi3{r,t) + df:iga{r,t) - {2/d)Sai3dxgx] - ■ySafsdxgx , (8) 

(Xa(r,<)) = [ie+p)/p]go.ir,t)-Kdc.T{r,t) , (9) 

where p(r, i) and r(r, t) are the local pressure and temperature fields, respectively, p is the equilibrium pressure and 
e is the equilibrium energy density; v and 7 are the kinematic shear and bulk viscosities, respectively, and k is the 
thermal conductivity, d is the spatial dimension, and cra/3 is the macroscopic viscous stress tensor. There are both 
non-derivative, reactive, and dissipative contributions to the constitutive relations. The form of the reactive terms 
can be inferred from Galilean invariance. The dissipative terms follow from an expansion of the current densities to 
first order in the gradients of the local conjugate forces g(r, t), p(r, t), and T(r, t); symmetry dictates the general form 
of these terms. Non-linear terms in the constitutive equations have been omitted because we only require the linear 
hydrodynamic equations in the following. 

The local equilibrium averages of Eqs. ©-(jSl), together with the constitutive relations given by Eqs. Q-® provide 
a complete description of the hydrodynamics of the liquid. The resulting linearized Navier-Stokes equation is 

dtga{r,t) + dap{r,t) - iydlga{r,t) - (^-f + dadxgx{r,t) ^ . (10) 

The corresponding equations for the mass density and energy density are 

dtpir,t)+d0f,{r,t)^O, (11) 

and 

5te(r,t) + [ie+p)/p]dxgx{r,t) - KdlT{r,t) = 0, (12) 

respectively. 

Because of the cell structure introduced in SRD to define the collision environment, angular momentum is not 
conserved in a SRD collision. As a consequence, the macroscopic viscous stress tensor is not a symmetric function of 
the derivatives dagp, and instead of Eq. (|SJ|, the constitutive equation has the general form |0 

TQ/3(r,t) =p(r,t)(5a/3 - vi[dagi3{r,t) + di3ga{r,t) - {2/d)Safidxgx] (13) 
- V2idi3gcir,t) ~ dcgp{Y,t)) - -fScpdxgx, (14) 

where V2 is a new viscous transport coefficient associated with the non-symmetric part of the stress tensor. Because 
the kinetic contribution to the microscopic stress tensor is symmetric, = and v^™ = It is also easy to 

show that 7*^™ = 0, so that the kinetic contribution to the macroscopic viscous stress tensor is 

al"f,\v,t) = ;.^-'"[9„gMr,i) +a^.9a(r,i) - {2/ d)5^pdxgx{vM (15) 

In Ref. [22 , it was also shown that the coUisional contributions to the viscous transport coefficients fulfill the relation 

[{d - 2)/d]<' - + 7'°' = 0, (16) 

and that the collision contribution to the macroscopic viscous stress tensor is 

<^'(r, t) = « + ^T')dfi9o.{v, t) = ;.^°'a^5a(r, t) , (17) 

up to a tensor with vanishing divergence, which will therefore not appear in the linearized hydrodynamic equations. 
The resulting linearized hydrodynamic equation for the momentum density is, therefore 

atg„(r,i) +a„p(r,t) - vdlg^{v,t) - ^j.'=™a„5A5A(r, i) = 0, (18) 

where v — ly^™ + v'^°^ and v^"'' and v'^"^ are the kinetic and collision contributions to the shear viscosity. The equations 
for the mass and energy densities remain unchanged. Comparison of Eq. H18|l with Eq. 110|l shows that the only 
difference between the Navier-Stokes equation for an isotropic liquid and an SRD fluid is in the coefficient of the 
dadxgx{^,t) term, where i/ is replaced by jy'^™. The bulk viscosity does not appear in Eq. I|18|l because it is zero 
for the SRD algorithm. Note that both equations are identical in d = 2 and that the only difference in d = 3 is a 
correction to the sound attenuation coefficient associated with the viscous dissipation of longitudinal sound waves. 
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IV. DYNAMIC CORRELATIONS 

Spontaneous thermal fluctuations of the density, p(r,t), momentum density, g(r, i), the energy density, e(j,t) are 
dynamically coupled, and an analysis of their dynamic correlation functions in the limit of small wave vectors and 
frequencies can be used to determine a fluid's transport coefficients. In particular, because it is easily measured in 
dynamic light scattering, x-ray, and neutron scattering experiments, the density-density correlation function — the 
dynamic structure factor — is one of the most widely used vehicles for probing the dynamic and transport properties 
of liquids 01 . 

In the following, we summarize the predictions of linearized hydrodynamics for the dynamic correlation functions 
of simple liquids in the hydrodynamic regime, and then use the results to analyze SRD simulation data in order 
to validate the theoretical results for the transport coefficients given in Table I. In particular, we provide in this 
way the first direct confirmation that the bulk viscosity is indeed zero for this model. Our discussion follows 
closely that of Ref. p7| . but is included because the derivation using Laplace transforms is not widely used in 
the literature, and the detailed results are needed in the subsequent analysis. The starting point is the linearized 
hydrodynamic equations given in Eqs. ^UJ), ((TT|l . and l(T^ . There are four modes, one transverse shear mode, and 
three coupled longitudinal modes. In order to keep the analysis general, we include the bulk viscosity in this section. 

Transverse fluctuations: Divide the momentum density g(r, t) into transverse and longitudinal components, 

g(r,<)=g||(r,i)+g^(r,i), (19) 

where V x g|| (r, t) = and V • g±(r, t) = 0. Taking the curl of the Navier-Stokes equation, the transverse component 
of the momentum density, gj^(r,t), is found to satisfy the diffusion equation 

dtg±{r,t) ^ iydlg^(r,t). (20) 

By performing the Fourier-Laplace transform (Im z > 0) 

g^(k,z) = rdte'^' f dv e-^'' '-gj.(r,t) , (21) 
Jo Jv 

the solution of the initial value problem, which describes the response of the transverse mode to an initial perturbation 
(5g_L (k, t — 0) from equilibrium, is 

r!g^^gi(k,z) <5gi(k,t = 0), (22) 

where = z + ik^v. 

Longitudinal fluctuations: For the longitudinal components, it is convenient to introduce the variable q{r,t), which is 
(T times) the entropy density, 

q{r,t)^e{r,t)-'-±Pp{r,t) (23) 
P 

in place of the energy density and use the relations 

Vp(r, t)=[^)^ Vp(r, t) + ^(^^^ Vg(r, t) (24) 

and 

. ^ fdT\ , , V fdT\ , , 
VT(M).(^-j^Vp(M) + -(^-j^V,(r,t), (25) 

where S is the total entropy, to eliminate the pressure and temperature fields. Taking the Fourier-Laplace transform, 
the resulting coupled set of equations for the longitudinal modes can be written as 
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where, for example, 5p(k, < = 0) = p(k.,t = 0) — p; p is the equiUbrium density, and k = |k|. Note that since 
(5g||(k, z) II k, the equation for the longitudinal component of the momentum density is a scalar equation. When 
writing Eq. (j26|l . we have used the relations 



T fdS\ 



P""^- V \dT) ' 



T fdS\ 



and 



dp 
dp 



dp 
dp 



(27) 



to simplify the final expression, c is the adiabatic speed of sound, and Di = 2[{d— l)/d]i' + 'j. Eq. H26(l describes how 
the longitudinal modes relax in response to initial perturbations 6p(k,t = 0), 6g\\(]i,t = 0), and dq(k,t = 0). The 
zeros of the determinant of the coefficient matrix, fJ, give the complex frequencies of the hydrodynamic modes of the 
system. For small wave vector k, the solutions of the resulting cubic equation are (up to terms of order k^) 



and 



z^±ck- -k^V 
2 



-ik^Dr 



(sound poles) 



(28) 



where Dt — K/{pCp) is the thermal diffusivity; T = Dt{cp/cv — 1) 
deriving Eqs. (|28|l and H29|) . the thermodynamic relation 



(heat pole), (29) 
Di, is the sound attenuation coefficient. In 



Drc^icp/c-u - 1) 



V f dT 



T \ dp 



dp 
dS 



(30) 



has been used. 

Correlation functions: The matrix of dynamic correlation functions 



is given by [22 



5,,(k,z)= / d{t~t') / d(r-r')e^^(*-*''-*-("-'-'^([A,(r,i)-(A,)][A,(r',0-(A,)]) 



S'y (k, z) = ifcsT(r2 ^)iixij{'k), 



(31) 



(32) 



where there is a sum over repeated indices. For the transverse modes, $7 is simply a scalar function defined in Eq. 
(|22|l : for the longitudinal modes, however, it is a matrix, and the subscript indices denote the modes Sp{'k,t = 0), 
6g\\ (k, t — 0), and (5(7(k, t = 0). XijO^) is the static susceptibility matrix. 

The correlation function for the transverse mode follows from Eqs. I|22|l and (|32|1 and Xg±g±0^) — P: ^^^d is given 

by 



g±9± 



{k,z) 



iksT p 
z + ik'^v 



Taking the inverse Laplace transform, 



Sg^g^{k,t) = pksTe 



-uk-'t 



For the longitudinal modes, inverting and using the results Xg||p(k) = 0, Xgygy (k) = p, limfc_>o Xpp(k) 
limfc„+o Xgq(k) = pTcp, and limfc_^o Xp<?(k) = {T /rn){dp/dT)p (where m is the particle mass), one finds 



Spp{k, z) = iksTp ( ^ 



T L 



z + ik^ {r + DT[cp/c^ - 1]) 
z^ — c^k^ + izk'^T 



Cp 



1 



z + ik'^Dr 



and 



Sqq(k, z) = iksT 



pCpT 



z + ik'^Dr 



(33) 

(34) 
p{dp/dp)T, 

(35) 
(36) 



for the two scalar modes. Note that all non-vanishing static susceptibilities are symmetric in fc, so that fc-dependent 
corrections to Xij 0{k^) and therefore negligible in these expressions. For details, the reader is referred to 
Refs. illsi. 
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The complete spectral transform of the time dependent density correlation function, the dynamic structure factor, 
Sppik, uj), is obtained by setting z ~ lo -\- i5 and taking the limit ^ — > in Eq. H35|) . The final result is [2^ . 



Spp{k,uj) = 2kBTp 



dp 
dp 



(C„/Cp)c2fc4r 



(1 - Ou/Cp)k^DT 



(t^2 _ c2fc2)^ + {ivk^ry Cj2 + (k^DT) 



c'k'] k^Dr 



(cj2 



{iuk^TY 



(37) 

In experiments, it is genenerally not possible to measure Spp{k,t) directly. In simulations, however, both Spp{k,t) 
and Spp{k,uj) can be measured for a range of mean free paths and collision angles. The simplest way to determine 
Spp{k,t) is to take the inverse Laplace transform of Eq. (|35|) . In light of Eqs. H28|l and H29|l . it is sufficient to keep 
terms 0{k) in real parts and 0{k^) in the imaginary parts when evaluating the resulting contour integral. The final 
result is 



Spp{k,t) = 2kBTp 



dp 
dp 



cos (ckt) 



^ - 1 ) Dt ) - sin (ckt) 



1-- le 

c 



Using Eq. H36|l . it can be shown that the correlation function for the entropy density, q{r,t), is given by 

Sggik,t)=pCpkBT^e-''^''"'. 



(38) 



(39) 



Note that Eq. H39|l provides an independent way to directly measure Dt- 

These results remain valid for the SRD fiuid. The only modification is that De = 2 [{d — l)/d] v^"'' + i^^"', so that 
the sound attenuation coefficient is 



Cv 



1 



1 



(40) 



Note that in two-dimensions, the sound attenuation coefficient for a SRD fiuid has the same functional dependence 



on Dt and v 



.kin 



,col 



as an isotropic fiuid with an ideal gas equation of state (for which 7 = 0). Finally, since 



SRD describes an ideal fluid, p = pkBT/m and Cp = kB/m ■ 



{d + 2)kB/2m. 



V. MEASUREMENTS 



In our SRD simulations in two dimensions, the mass, momentum, and energy densities are measured at the cell 
level. The cell densities, A'^(^, i), are defined at the discrete set of coordinates ^ = am, with mp = 1, . . . , L, for each 
spatial dimension TT]. A superscript c will be used to denote that the corresponding quantity is defined at the cell 
level. The Fourier transform of the cell variables are 

A^(k,i)=^A^(|,t)e'''-S (41) 
s 

and the inverse transform is 

k 

The Fourier-Laplace transforms of the corresponding dynamic correlation functions are 

/>oo 

Sf^iKz)^ d{t^t')Y,e^<'~''^~^^-H[AUtt)~{A'l)W,{0,t')-{A'^,^ . (43) 
^0 J 

Transverse fluctuations: Instead of the evaluating the correlation function of the transverse component of the momen- 
tum density, it is more convenient in simulations to measure the vorticity, w(r,<) = V x g^(r,t). In two dimensions 
the vorticity is a scalar, Wz(r, t), and the dynamic correlation function decays as 

5^^^^ (fc, t) = eS^^g^ (k, t) = pkBTk^e--^'K (44) 

Simulation results for the normalized vorticity correlation function for A/a = t ^/kj^Tjrn — 1.0 with collision angle 
a — 120° and A/a = 0.1 with a = 60° are shown in Figure ^ Here, as with all results presented in this paper, 
averages are taken over 400,000 iterations and five different random number seeds. The solid lines in Figure ^ are a 
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plot of Eq. H44() using the theoretical prediction for the shear viscosity (sum of kinetic and collisional contributions) 
given in Tabled for the smallest wave vector k = (27r/L)(l, 0). The agreement is excellent. We have also fitted the 
decay profiles for the lowest two wave vectors, namely to k = (27r/L)(l,0) and k = (27r/-L)(0, 1) with Eq. I|44|l . and 
averaged the result to obtain estimates of the shear viscosity as a function of collision angle a for different mean free 
paths. The results are presented in Figure |2 and, as expected, the agreement between measured viscosities and the 
expressions given in Tabled is very good. These measurements clearly show that the theoretical expressions for the 
shear viscosity are accurate even for intermediate mean free paths. 

Longitudinal fluctuations: Density fluctuations were measured at the cell level, and their Fourier-Laplace transform 
is taken to determine the structure factor. A naive implementation of this procedure gives wrong results in the large 
frequency region of the spectrum J4J ■ resulting in finite contributions at all frequencies. This problem is well known to 
experimentalists [32| , the solution is to first do a Fourier transform to obtain density-density correlations as a function 
of time, symmetrize this result around t = 0, and then perform the Laplace transform from — oo to oo. Results for 
the structure factor for A/a ~ 1.0 with collision angle a = 120° and A/a = 0.1 with a = 60° are shown in Figures 
ISfa) andl^Jb), respectively. The solid lines are the theoretical expression given by Eq. H37|) using c = y/ dkBT/m 
and values for the transport coefficients obtained using the expressions in Tabled assuming that the bulk viscosity 
7 0. The agreement is excellent. There are three Lorentzian peaks, a central "Rayleigh peak" caused by the heat 
diffusion and two symmetrically displaced "Brillouin peaks" caused by the sound waves. The dotted vertical lines 
in the figures show the theoretically predicted frequencies of the adiabatic sound waves in a fluid with an ideal gas 
equation of state. 

We have also measured time-dependent density correlations for various wave vectors. Figures0fa) and^fb) contain 
a comparison of the measured time-dependent density correlation functions with the predictions of Eq. H38|) . for 
A/a = 1.0 with collision angle a = 120° and A/a = 0.5 with a = 90°. The agreement is excellent for all wave 
vectors considered. The bulk viscosity and thermal diffusivity were also independently measured by fitting these time 
dependent density correlations to the form given by Eq. H38|) while using the theoretically predicted shear viscosity 
in the sound attenuation coefficient, keeping Dt and 7 as free parameters. The results are shown in Figure [S] as 
a function of the wave vector squared, for the same set of parameters as in Figure 0] Once again, the theoretical 
expression for Dt is confirmed, and the bulk viscosity is indeed zero. 

In order to obtain an independent measure of the thermal diffusivity, we have measured the temporal behavior of 
the entropy correlations, S^^{k,t). The results are shown in FigureElfor A/a = 1.0 with collision angle a = 120° and 
A/a — 0.5 with a = 90°. As expected, these corrrelations decay exponentially for all wave vectors considered. The solid 
lines in Figure|Hlare a plot of Eq. (|39|l for the smallest wave vector k = (27r/L)(l, 0), using the theoretical prediction 
for the total thermal diffusivity, Dt (see Tabled, and the agreement is again very good. As was done for the vorticity 
measurements, we have also fitted the decay profiles for the lowest two wave vectors with Eq. (15^ . and averaged the 
results to obtain independent measurements of the thermal diffusivity. The results for these measurements are shown 
in Figure[7|as a function of the collision angle a for A/a = 0.5 and A/a = 1.0. The theoretical values obtained using 
the formulae for Dt in Table (sum of kinetic and collisional contributions) are shown in solid lines. These results 
are the first direct equilibrium measurements of the thermal diffusivity. 

Finally, it is important to emphasize that just as for the shear viscosity, collisional contributions provide the 
dominant contribution to the thermal diffusivity at small mean free path. Figure |S1 shows the theoretically predictions 
for both the collisional and kinetic contributions to the shear viscosity and thermal diffusitivity (inset) as a function 
of the mean free path A/a. Collisional contributions to both transport coefficients are particularly important for small 
mean free paths and small M . 



VI. CORRELATION EFFECTS 

Green-Kubo relations for the SRD transport coefficients have been derived in Ref. ^3 ^^'^ analyzed in Refs. 
[SIM El 113, where it was shown that there are both kinetic and collisional contributions to the shear viscosity 
and the thermal diffusivity. The collisional contributions to these transport coefficients have been discussed in detail 
in Ref. 

The kinetic contribution to the transport coefficients have been derived by several groups P, IE US 13 assuming 
molecular chaos. The results of these calculations are summarized in TableQ] Simulation results for the shear viscosity 
and thermal diffusivity have generally been found to be in good agreement with these predictions. However, it is 
known that there are correlation effects for A/a smaller than one They arise from correlated collisions between 
particles that are in the same collision cell for more than one time step. In the following, we expand on the discussion of 
Ref. and calculate the first correlation corrections to both the viscosity and the self-diffusion coefficient explicitly. 
Similar calculations can in principle be done for the correlation contributions to thermal diffusivity. The reason that 
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these corrections to v and Dt are generally negligible is that they are only significant in the small A/a regime, where 
the collisional contribution to the transport coeiEcients dominates. 

Figure IHl shows a comparison of both kinetic and collisional contributions to shear viscosity and thermal diffusivity 
(inset) for M — 3 and a = 60°. Correlation efi^ects would be most pronounced when both contributions are comparable, 
i.e. at a mean free path of A/a ~ 0.25 for v and \/a •2^ 0.1 for Dt (see Figure|H|. On the other hand, because there 
are no collisional contributions to the self-diffusion coefficient, correlation corrections dramatically increase the value 
of this transport coefficient in the small A/a regime. It is important to note that there are no correlation corrections 
to the collisional contributions so that the expressions for the collisional shear viscosity 5, 21, 22] and collisional 
thermal diffusivity IIJ,!!^! are exact. 

In the following, we restrict ourselves to two dimensions; the same analysis, however, can also be used in three 
dimensions. Expressions for the shear viscosity and the self-diffusion coefficient in this section obtained in the molecular 
chaos approximation will include contributions from fluctuations in the number of particles per cell. However, when 
calculating correlation corrections, we will assume that the number of particles per cell, M, is fixed. Including these 
fluctuations is straightforward but tedious, and since it would not provide any additional insight into the underlying 
phenomena, we have decided to ignore this effect in the following. 



A. Shear viscosity 

The Green-Kubo relation for the kinetic contribution to the shear viscosity is 

oc 

^i^^n ^kBTT^'Ginr) (45) 

n=0 

where the prime indicates that the t — Q contribution to the sum occurs with the weight 1/2, and 

G{nT) = ^(w„(0)z;,2,(0)wj,(nT)w,j,(nT)}/Ar(fcsT)2 . (46) 

In the molecular chaos approximation j2ll |. 

G(nr) = Gcinr) = [l - 2 sin^ {a) {M - 1 + e"*^) /M] " . (47) 

Inserting this expression into Eq. I|45|l and summing, one obtains the kinetic contribution to viscosity given in Table 
n Figure 13 contains a comparison of simulation results for Ginr) with the molecular chaos approximation Eq. I|47|l . 
As can be seen, the flrst and the most important correlation contribution to v^"'' occurs for n = 2. The functional 
form of this leading correlation correction, 5G{2t) = G{2t) — Gc{2t), can be calculated analytically. 

As illustrated schematically in Figure [TUI there are six distinct particle configurations which contribute. The first 
two, shown in Figures lior i'l and I10r 2') occur when i = j in the sum in Eq. H46() : we will call these the diagonal 
contribution. In Figure ITUT l). particle k is in the same collision cell as i for both t = 2t and t = t. In Figure [TUT 2). 
two distinct particles, labeled k and I, are in the same collision cell as i at both t = 2t and t = t. Other (off-diagonal) 
contributions, which occur for i ^ j, are given in Figure^j3)-(6). These contributions are significant only at small 
mean free paths, since their amplitudes are proportional to the probability that two or more particles are in the same 
collision cell for multiple times. 



1. Diagonal contributions 

The first diagonal contribution, which we denote by 6Gi(2t), occurs when two particles, with indices i and k, are 
in the same collision cell at both t — t and t — 2t. The probability for this to occur is p2; P2 is calculated in the 
A/a — > limit in Appendix A. In two dimensions the velocity of particle i at time t — nr is related to the velocities 
of its collision partners at t = (n — l)r by 



Vixinr) = cvix{[n - 1]t) + ^^'^Vkxi[n - 1]t) + s Viy{[n - 1]t) - —'^Vkyi[n - 1]t) 



(48) 



Viyinr) = cviy{[n ~ 1]t) + 



l__c 

M 



^Vky{[n - 1]t) 



w„([n - l]r) - j-^^Vkx{[n - 1] 
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where c = cos(a) and s = siii(a). Using Eq. (|48|l to relate the velocities at t = 2t to those at < = r, we have 

6G,{2t) = 2p2{M - l)Y,(v^40)^^viQ) [Civ:4r) + C2V^y{r)] (^^Vky{r) + ^t^.^r)^ ^ /7V(fcBT)2, (49) 

where Ci = 1/M + c(l — 1/M) and (2 = s{l — 1/M), and the factor (M — 1) accounts for the sum over k i and 
the factor 2 comes from the fact that i and k can interchange roles. The equilibrium average in Eq. (|49|) entails the 
average over all initial coordinates and velocities (at t = 0) as well as averages over the stochastic rotations (±a) at 
t = T and t = 2t. Performing the average over the collision angle at t — 2t in Eq. (|49|l removes all terms linear in s, 
so that 

SGi{2t) = ^(^;„(0)7;,,(0)[Ci(l - c)v,^{T)vky{T) + C2Sv^y{T)vkAr)]) /NikeTf . (50) 

i 

Using Eq. (|48|1 once again with n — \^ 



5G,{2t) = ?M^-^^(^;„(0)t;.,(0) 



Cl(l - C)[C1V,.{0) + C2V^ym (^^^yi^) + ]^^«(0) 

+ C2s[Ci«^y(o) - C2f«(o)] (^^-(0) - 



)/N{kBTf. (51) 



Averaging now over the collision angle at t = t and the particle velocities and coordinates at i = yields 

5G^{2r) = ^^^^^ ~ [Ci(l - c) + Q^sf . (52) 

The other diagonal contribution, which we will denote as 5G2{2t), arises when three particles, with indices i, k and 
I, are in the same collision cell at t = t and t = 2t. The probability that three particles are in the same collision cell 
in consecutive time steps will be denoted by ps; ps is calculated in the A/a ^ hmit in Appendix B. Using Eq. I^HJ, 

SG2{2T)=p3{M~l){M-2)Y,(^v.40)v^y{0) (l^VkAr) - {^Vky{r)^ (^^viy{r) + ^viAr^j'^ /NiksTf , 

(53) 

where the prefactor (M — t){M — 2) comes from the double sum over k, I. Averaging over the collision angle at i = 2r 
yields 

SG2{2r) = P3(M-^KM-2) ^ ^^^^(o)„^^(o) ([^ _ c]^v,Ar)v,y{r) - s\,y{r)vUr))) /NiksTf. (54) 

i 

Finally, using Eq. (|48|l again and averaging over the collision angle at t — t and velocities and coordinates at i = 
yields 

P3(M-l)(M-2) ^ 2 
dG2(2T) = — [2c(c - 1)] . (55) 



2. Off-diagonal contributions 

The analysis of these contributions is very similar to that which was used to evaluate the diagonal contributions, 
so we provide fewer details than in the previous subsection. There arc four off-diagonal contributions, all of which 
contribute with probability p2. The first, (5G'3(2r), shown in Figure [TUT S). can be written as 



2 



<5G3(2r) = 2p2{M - l)(Af - 2) ^ /^v,AO)v^y{0) (^«- W " J^^^vi^)) (^^^2'^^) + ]^"^-(^)) ) /^C'bT) 

(56) 

where the factor 2(Af — l)(Af — 2) accounts for i and k interchanging roles and the double sum over j and k. Following 
the same procedure as was used to evaluate 6Gi{2t) and (56*2 (2t), we have 

SG,i2r) = 2p.(M-^l)(Af-2) ^^^^^ _ ^^^^^^^^ _ ^ ^^^^ ^^^^ 
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The second contribution, SG4{2t), depicted Figure 11^ 4) can be written as 
6G42r) = 2p2(M - 1)(M - 2) ^ ^t^..(0)«.,(0)[Ci^',.(r) + C2V,y{r)] (^i^^,,,(r) + ^VkAr)^ ) /N{kBT)^ (58) 

where the prefactor is similar to that of SG3{2t). Using this result, it is straightforward to show that 

SG.(2r) ^ ^"''"-jf^-^' lM. - DllCd - c) + . (59) 

The third contribution involves the configuration shown in Figure [TUT S). It is similar to the diagonal contribution 
(5G2(2r) except for the fact that the particle with index i is not in the same shifted cell as k and I at time t = 2t. 
Only two particles are in the same collision cell for consecutive time steps, so that the relevant probability is p2- By 
analogy with the expression for <5G2(2r), 

(0 \ P2(M~l)(M-2)(Af-3) ^^ 2 
(5G5(2t) = — [2c(c-l)] , (60) 

where the factor of [M — 3) comes from the additional sum over j with the constraints i ^ j ^ k ^ I. 

The final contribution, 5Gq{2t)^ occurs when i and j are in the same collision cell at both t — 2t and t — t (see 
Figure Ej6)). This contribution is given by 

<5G6(2r) = 2p2{M - 1) ^ ^i,„(0)«.,(0) [l^v^^r) ~ ^^^.y(r)) [Civ.yir) - C2t',.(r)]^ iNiksTf, (61) 

where the factor 2(M — 1) accounts for the sum over j and interchanging i and j. Following the same procedure as 
for the other contributions, we find 

5G,{2t) ?^al^l^[Ci(l - c) + . (62) 
The total correlation enhancement is obtained by summing these six contributions, 

6 

5G{2t) =^5Gn{2T) . (63) 

Note that the only dependence on the temperature and time step r occur in the probabilities p2 and p^. The measured 
correlation contributions to the shear viscosity are shown in Figure ^2 as a function of A/a for a = 60°, 90°, and 120°. 
Using the A/a ^ values for the probabilities p2 and ps calculated in Appendices A and B, one finds 

<5G(2r) ^ ^^"(2r) = [{4M[1 + cos(a)]}2 + 7(2 - M) cos'{a)] . (64) 

n— 1 



and 



,ki' 



V 



keTT {y"'Gc{nT) + 6G{2T)) (65) 



, Ji=0 



for the shear viscosity. As can be seen from Figure ^2 the results obtained using Eq. (|64|l (asterisks) are in 
excellent agreement with simulation data in the limit of zero mean free path. More generally, we have determined 
the probabilities p2 and p^ numerically and have found that they depend only on the value of the mean free path. A, 
and not on r and T individually (see inset of Figure IT^ for a plot of p2 as a function of A/a). Finally, these results 
can be used in Eq. H63|l to obtain an estimate of the correlation contribution 5G{2t) for arbitrary A. The asterisks in 
Figure El show the results of this procedure, and as might be expected, the agreement is excellent. Finally, simulation 
results for the total kinetic contribution to the viscosity as a function of time are shown in the inset to Figure El 
The filled squares (■) are the predictions of molecular chaos approximation, and asterisks are a plot of Eq. I|65|l . 
The incipient long-time tail is clearly visible in the figure. This is one of the reasons that it has been difficult to 
obtain good estimates for the "bare" kinetic contribution to the transport coefficients. In principle, the methods used 
to determine 5G{2t) can also be employed to determine these correlation contributions at greater time lags. The 
corresponding probabilities, however, that particles become collision neighbors after a finite time interval are much 
harder to determine, since they depend in detail to relative probabilities of various fluctuating flow configurations. 
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B. Self-diffusion coefficient 



The Green-Kubo relation for the self-diffusion coefficient is [3 



1 rri 



(66) 



n=0 



and 



H{nT) 



{v,{Q)-v,{nT))/{kBT) . 



(67) 



Since the self-diffusion coefficient is a single-particle property, there is no sum over i in Eq. H66|l . In the molecular 
chaos approximation, 



Inserting Eq. (|68|l into Eq. (|66|l and summing the resulting geometric series, one can obtain the expression for the 
self-diffusion coefficient in two dimensions given in Table n 

The self-diffusion coefficient is unique in that there is no "coUisional" contribution; as a result, correlation corrections 
are much more important at small mean free path and can lead to large corrections to the results obtained using the 
molecular chaos approximation. Corrections to this result occur when two or more particles occupy the same collision 
cell at different time steps. Figure El contains a comparison of the molecular chaos approximation for the velocity 
auto-correlation function, H{nT), (■) with simulation results (o). 

The first of these correlation corrections, 5H{2t), occurs aX t = 2t. The contributing configuration, in which two 
particles, i and fc, are in the same (shifted) cell at both t — t and t = 2r, is shown in Figure El the probability for 
this to occur is again p2 ■ The contribution of this configuration to the velocity auto-correlation function is 



where the factor 2 arises since both x and y components contribute; the factor M — 1 accounts for the fact that k i 
can correspond to any of the M — 1 particles. Following the procedure outlined in the discussion of correlation effects 
to the viscosity, it is straightforward to evaluate the averages in Eq. (|69l) . The final result is 



This is the only correlation correction at i = 2r. At longer times there are similar higher order correlation effects 
arising, for example, when two particles are in the same shifted cell for three time steps, etc. It is straightforward 
but tedious to calculate these contributions. The probability p2 in Eq. I|7U|) is determined analytically in the limit 
A/a — > in Appendix A, where it is shown that p2 — (2/3)'' in d dimensions. For finite A/a, it can be measured 
in simulations. It should be noted that p2 is related to the quantity Ci of Ref. ^26J, which denotes the number of 
particles that are neighbors of a given particle for two consecutive time steps; more precisely, p2 = (Ci ~ l)/(-^^ ^ !)• 

Figure is a plot of simulation results for SH{2t) as a function of A/a for three different values of the collision 
angle a. The asterisks (*) are result of Eq. (|70|l using the A/a — + prediction p2 = 4/9; as can be seen, the agreement 
with simulation data is excellent. The inset in Figure shows p2 as a function of A/a. The value for the probability 
P2 in the limit A/a —> is in excellent agreement with the result derived in Appendix A. 

We have only considered correlation effects caused by two particles occupying the same collision cell for the two 
consecutive time steps. In fact, additional contributions arise any time two particles find themselves in the same 
collision cell for more than two time steps or after any number of time steps. It is interesting to see just how 
important these latter contributions are by summing up all possible contributions from 5H{2t) and Hc{t). Figure 
II 51 shows the contributions for the first five time steps in the series. In this approximation, the self-diffusion constant 
can be written as 




(68) 




(69) 




(70) 




(71) 



where 



5H{nT) 



0, 

5H{2t), 
5H{2t)H,{t), 
{l/2)5H{2Tf + {3/2)SH{2t)Hc{2t), 
{3/A)SH{2t)'^Hc{t) + 2SH{2t)Hc{3t) 



n<l 

n = 2 

n ~ 3 

n = 4 

n — 5 



(72) 
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are the contributions shown in Figure^] The solid Unes in the figure represent the factor SH{2t), and the dashed fines 
represent the factor Hc(t) from time steps where particles are not correlated and the molecular chaos approximation 
is valid. The resulting contribution of Eq. 1)72(1 is shown in FigureElby asterisks (*). The agreement with simulations 
at i = 2t is perfect, as expected, and the prediction for larger time intervals is improved. It would be extremely 
interesting to extend the analysis to consider the contributions of correlation affects at larger time intervals that lead 
to long time tails in the velocity auto-correlation function [s^ . 



This paper contains the first detailed analysis of equilibrium dynamic correlation functions using the SRD algorithm. 
The dynamic structure factor, vorticity and entropy density correlation functions were measured and used to provide 
unbiased estimates for the viscosity, thermal diffusivity, and bulk viscosity. The results are in good agreement with 
earlier numerical and theoretical results, and provide the first direct verification that the bulk viscosity is zero for this 
algorithm. 

Table ^ contains a complete summary of analytical results for the transport coefficients of this model, and the 
results of this paper verify that we have an excellent understanding of the SRD algorithm at the kinetic level and that 
the analytic expressions for the transport coefficients do indeed provide a very accurate description of the SRD fluid. 
Furthermore, the analysis of the dynamical structure factor and the dynamic entropy density correlation function 
verify directly that the algorithm satisfies the fiuctuation-dissipation theorem. Verification studies of this type will 
be important for generalizations of the algorithm which model excluded volume effects through the use of biased 
multi-particle collision rules which depend on the local velocities and densities. 

Finally, corrections to the self-diffusion coefficient and shear viscosity arising from the breakdown of the molecular 
chaos approximation at small mean free paths were analyzed. In addition to deriving the form of the leading correlation 
corrections to these transport coefficients, the probabilities that two and three particles remain collision partners for 
consecutive time steps are derived analytically in the limit of small mean free path. Extensions of this approach could 
be used to study the development of long time tails in the velocity and stress auto-correlation functions. 



Support from the National Science Foundation under grant No. DMR-0513393 and ND EPSCoR through NSF 
grant EPS-0132289, is greatfuUy acknowledged. 



Random grid shifts in x- and j/-directions are statistically independent. We therefore first calculate the contribution 
from shifts in the x-direction; the final probability for general shifts in two-dimensions is then obtained by squaring 
this one-dimensional result. The following calculations are done in the limit A/a ^ 0, so that the particles do not 
move between time steps. 

There are three different ways for two particles to be in the same shifted cell at consecutive time steps (see Figure 
I16|l . If the shifted cell index at i = t, £,s{t) is zero, then ^s(2t) can be either —a, or a. Because these are mutually 
exclusive events, the probability of each scenario has to be calculated separately and then summed. 



VII. CONCLUSION 
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APPENDIX A: CALCULATION OF pa 



1- C.(^) = 0, Cs(2r)=0 



The situation is shown pictorially in Figure [TBT al 
consecutive time steps can be written as 



The probability P2 



<2 that two particles are in cell = for two 




X 



e {X2 - 62) [l-e{x2-a- 62)] , 



(Al) 
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where 5i and S2 are the shifts at times r and 2t, respectively. Making the substitutions Xi = xi—Si and X2 = X2—S1, 
Eq. (|A1I) becomes 



— d62 dSi / dXi / e {Xi + Si - S2)[l - e {Xi - a + Si - S2)] 

^ J-a/2 J-a/2 JO JO 

X e(X2 + <5i-<52)[l-e(X2-a + 5i-(52)]. (A2) 



To simphfy further, introduce pi = 61— S2 and integrate Xi and X2 over the portion of the square where the integrand 
is non-zero. This yields 



P2^^ 1 dS2 I dSi{a-\pi\f . (A3) 



a/2 i-a/2 
dS2 

-a/2 J-a/2 



Using, finally. 



Eq. (|A3|I gives 



ra/2 /■a/2 

dS2 / dSi = a^, (A4) 

-a/2 J~a/2 

/.a/2 i.a/2 3 

/ d52 / d^llpil = - (A5) 
J -a/2 J-a/2 

/•a/2 /-a/Z 4 

/ dS2 / rf<5ip? = -, (A6) 

J-a/2 J-a/2 » 



= !- - + - = -• (A7) 
^2 3 6 2 ^ ' 



2. Cs(r) =0, C.(2r) = -a 

This situation is illustrated in Figure ITHb. The probability P2 that two particles are in cell = at i = r and 

= —a at t = 2r is 



1 ra/ ^ r^i ^ ra^oi /•a-toi 

P2=— dS2 dSi dxi dx2 e{xi+a-S2)[l-e{xi-52)] 

O- J-a/2 J-a/2 J Si J Si 

X e(x2+a-(52)[l-e(x2-(52)] . (A8) 
Making the substitutions Xi = xi — Si and X2 = X2 — Si, Eq. (jA8|) becomes 

ra/2 pa/2 

dS2 / 

2/2 J-a/2 



pi 



pa/2 pO'/Z pa pa 

— dS2 dSi / dXi / dX2 e{Xi+a + Si~S2)[l~e{Xi+Si~S2)] 

^ J-a/2 J-a/2 JO JO 

X Q{X2 + a + Si-52)[l-Q {X2 + <5i - 62)] . (A9) 



As in the previous subsection, introducing pi = Si — S2 and integrating Xi and X2 over the portion of the square 
where the integrand is non-zero yields 

Pf = ^ f'^ dS2 rf<5i|pipe (-Pi) = . (AlO) 

a* J-a/2 J-a/2 12 

3. 5s(t) = 0, ^.(2r) = a 

Referring to Figure 116b . the probability that two particles are in cell = at i = r and = a at i = 2t is 

2^ l'a/2 p(2/2 pa+Si pa+Si 

P2 = ^ dS2 dSi / dxi / dx2 & {xi - a + S2) [1 - <d {xi - 2a + S2)] 

J-a/2 J-a/2 J 61 J Si 

X Qix2-a + S2)[l-e ix2 -2a + S2)] . (All) 
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Making the same change of variables as in the previous two cases, Eq. IjAlip becomes 

P2^— dS2 dSi dXi dX2 e{Xi-a + 6i + d2)[l-e{Xi-2a + 5i+S2)] 

^ ■/-a/2 J-a/2 JO JO 

X e{X2-a + di + S2)[l~eiX2-2a + Si+S2)] . (A12) 
Introducing now p2 = Si + S2 and performing the integrals over Xi and X2 yields 

^ ' dS2 — " — 

a/2 J~a/2 



P2=- dS2 dSWQ (P2) = -T . (A13) 

./_„/9 ./_„/9 12 



The final result in two dimensions is obtained by summing the results given in Eqs. HA7|I . HA10() and HA13|I . and 
squaring, so that 

P2 = {p^+V^+V^f = \- (A14) 

APPENDIX B: CALCULATION OF pa 

The calculation of ps in the limit A/a — > is similar to that of p2 in the previous Appendix. There are three 
scenarios, as depicted in FigureEI(with three particles instead of two), is the probability that three particles are 

in the = for two consecutive time steps: 



= — I d52 I d6i I dxi / dx 



a/2 pa/ 2 pa-\-6i pa-\-5i pa-\-Si 

C2 



dx3 8 {xi 


- '^2) [1 


— Q {xi — a 


-S2)] 


X e {x2 


- 62) [1 


— Q {x2 — a 


-S2)] 


X 6(2:3 


- S2) [1 


— Q {x3 — a 


~S2)] 



(Bl) 

To evaluate this integral, make the same change of variables to Xi and X2 as in the previous Appendix and introduce 
Pi — Si ~ 62- Performing the X integrals then gives 

p^ = \ [ ' dS2 ( ' dSi (a - IpiD' . (B2) 



a" 



a/2 J-a/2 



Using 

fa/2 na/2 

dS2 I d<5i|pi|3 = r_ (B3) 

! 

and Eqs. (HU, (HsJ, ^M^ and Eq. ^ yields 



-a/2 J-a/2 10 



pt = l-l + -- — = -. (B4) 
2 10 5 ^ ^ 

The calculations of pf and p§ are similar to those outlined in Sections lA 21 and lA 31 and both are equal to 1/20. 
Summing these results, 

P3^{pi+Pi +P^f = \ (B5) 

in two dimensions. 
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TABLE I: Theoretical expressions for the shear viscosity u, the thermal diffusivity Dt, and the self-diffusion coefficient D, in 
both two and three dimensions. M denotes the average number of particles per cell, a is the collision angle, fc^ is Boltzman's 
constant, T is the temperature and r is the time step. Except for self-diffusion constant, where there is no coUisional contribu- 
tion, both the kinetic and coUisional contributions are listed. The expressions for shear viscosity and self-diffusion coefficient 
include fluctuation corrections for small M; however, for brevity the relations for thermal diffusivity are correct only up to 
0(1/M) and 0{\/M^) for the kinetic and coUisional contributions, respectively. Both kinetic and coUisional contributions to 
shear viscosity have been calculated using two complementary approaches, equilibrium Green- Kubo relations I19II20II21L 
and a non-equilibrium approach 5, 18]. Results from both approaches are in complete agreement. Similarly, the kinetic con- 
tribution to thermal diffusivity has also been calculated using these two approaches. The predictions of both approaches are 
identical in two dimensions and agree up to (and including) 0{1/M) in three dimensions; higher order contributions in 1/M 
were not considered in the Green-Kubo approach. The kinetic contribution to thermal diffusivity calculated using the non- 
equilibrium approach was taken from Ref. since there appear to be several misprints in Eqs. (62), (63) and (64) in Ref. 
[3. The coUisional contribution to thermal diffusivity has been calculated in both two and three dimensions in Refs. |2ll |22| . 
Because of space limitations, only the leading terms in 1 /M are given here. For the complete expression, the reader is referred 
to |23 |. To the best of our knowledge, the fluctuation corrections for self-diffusion coefficient are presented here for the first 
time. 



Transport coefficient 


Dimension (d) 


Kinetic (xfcsrr/2) 


CoUisional (xa^/r) 


Shear viscosity, p 


2 
3 


M 1 

(M-l-Fe-JW)sin^{ci) 

5 A/ 1 

(AJ-l + e-*-' )[2-cos(a,)-cos(2a)l 


^(M-l + e-«) [l-cos(a)] 


Thermal diffusivity, Dt 


2 
3 




3{<i+2)Af (1 if) [1 cos(a)] 


Diffusion coefficient, D 


2 
3 


dM 1 

(l-cos{a)){Af-l + e-J« ) 
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FIG. 1: Normalized vorticity correlations as a function of time for k = ^(1,0) (dotted-dashed lines), k = ^(1, 1) (dotted 
lines) and k = %-(2, 0) (dashed lines). The solid line shows Eq. 14411 using the theoretical expressions given in TableQlfor shear 
viscosity. The decay profiles were fitted to Eq. 1441 to obtain values for the shear viscosity u. a) A/a = 1.0 , a = 120°, b) 
A/a = 0.1, Q = 60°. Parameters: L/a = 32, M = 15 and r = 1.0. 
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FIG. 2: Shear viscosity z/ as a function of collision angle a for mean free paths A/a = 0.1 (►),0.5 (■), and 1.0 (•) measured 
using the decay of the vorticity correlations. The solid lines are the theoretical prediction given in Table El i.e. the sum of the 
kinetic and coUisional contributions. Parameters; L/a = 32, M = 15 and r = 1.0. 
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FIG. 3: Normalized dynamic structure factor, Spp{k,uj)/xpp{k), for k — ^(1,1) and a) A/a = 1.0 with a — 120°, and b) 
A/a — 0.1 with a = 60°. The solid line is the theoretical prediction obtained using Eq. 13711 and the expressions for the 
transport coefficients given in Table |l] The dotted lines show the predicted positions of the Brillouin peaks using the dispersion 
relation u — ±cfc. Parameters: L/a = 32, M = 15 and t — 1.0. 
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FIG. 4: Normalized density correlations, Spp{k,t)/ Spp{k,0), as a function of time, for k = ^(1,0) (•), k = ^{1, 1) (■), and 
k = ^(2,0) (►). The solid lines are the theoretical predictions obtained using Eq. I38II and the expressions of transport 
coefficients given in Table|I] a) A/a = 1.0, a = 120°, and b) A/a = 0.5, a — 90°. Parameters: L/a = 32, M = 15 and r = 1.0. 
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FIG. 5: Fitted values for the bulk viscosity 7 (•) and thermal diffusivity Dt (□) as a function of the wave vector squared for a) 
A/a = 1.0 and a — 120°, and b) A/a = 0.5 and a = 90°. Data was obtained by fitting the time dependent density correlations 
(see Figure 3J using Eq. 1381 while keeping Dt and 7 as free parameters. Dashed and dotted lines represent the theoretically 
predicted values for 7 and Dt, respectively. System size L/a ranges from 32 to 128. Parameters: Al = 15 and r = 1.0. 
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FIG. 6: Normalized entropy correlations, Sqq{k,t)/ Sqq{k,0), as a function of time for k = ^(1,0) (dotted-dashed lines), 
k = ^(1, 1) (dotted lines), and k = (dashed lines). The solid line shows Eq. (I39II using the theoretical expressions for 

the thermal diffusivity, Dt, given in Tabled The decay profiles are also fitted to Eq. I|39|l to obtain unbiased estimates for the 
thermal diffusivity. a) A/a = 1.0 , a = 120°, b) A/a = 0.5, a = 90°. Parameters: L/a = 32, M = 15 and r = 1.0. 
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FIG. 7: Estimates for the thermal diffusivity Dt as a function of coUision angle a for mean free paths \/a — 0.5 (■) and 
1.0 (•) obtained by fitting the decay of entropy correlations, Sqq(k,t). The solid lines are the theoretical prediction given in 
Table 111 i.e. the sum of kinetic and coUisional contributions. Parameters: L/a — 32, M = 15 and r = 1.0. 
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FIG. 8; Shear viscosity v and thermal difFusivity Dt (inset) as a function of A/a. Both plots are obtained using the theoretical 
expressions given in Table |I] The solid and dotted lines are the kinetic and coUisional contributions, respectively. The dashed 
lines are the total contributions to these transport coefficients. For consistency, in the calculation of thermal diffusivity, both 
the kinetic and coUisional contributions are taken only up to and including 0(1/A/). In the plots. A/a was varied by changing 
ksT for fixed r = 1.0. Parameters: M = ?, and a = 60°. 
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FIG. 9: Stress correlations, G(nT), as a function of time step n for A/a = 0.25 and AI — 3. The inset shows the normalized 
kinetic contribution to shear viscosity, v*"^^ /ksTT, as a function of time step. The measured values are open circles (o), and 
the results of molecular chaos approximationi, Gcinr) (see Eq. I|47^ 'l. are filled squares (■). The asterisk (*) in the main figure 
is the prediction of Eq. I63II using the numerically determined values of the probabilities. The asterisks in the inset are a plot 
of Eq. jnSJ. Parameters: L/a = 64, a = 60°, r = 1.0. 
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FIG. 10: Diagrams contributing to correlations at t = 2t in the calculation of kinetic contributions to shear viscosity. The first 
two diagrams show the diagonal and the others show the off-diagonal contributions. 



28 




FIG. 11: Correlation contributions SG{2t) as a function of A/a. Results for collision angles a = 60° (o), 90° (□), and 120° (>) 
are presented. The asterisks (*) are the theoretical predictions in the limit A/a —^ 0. Parameters: L/a = 64, M = 5, r = 1.0. 
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FIG. 12: Normalized velocity auto-correlation function H{nr) as a function of time step n. Measurement values are shown by 
open circles (o), the geometric series by filled squares (■), and the sum of the geometric series and the correlation contributions 
fEa.l72l using the numerically determined P2) are shown by asterisks (*). The inset shows the normalized diffusion coefficient, 
D/ksTT, as a function of the time step n. The asterisks in the inset are a plot of Eq. (|nj 

using Eq. 17211 . Parameters: 

A/a = 0.1, a = 90°, L/a = 64, M = 5 and r = 1.0. 
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FIG. 13: Schematic diagram showing the configuration contributing to correlations at f = 2r in the calculation of self-diffusion 
coefficient. Particles i and k are in the same shifted collision cell at both t = t and t = 2t. There are M — 1 such contributions 
with probability p2- 
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FIG. 14: Correlation contributions at f = 2r as a function of A/a for the self-diffusion coefficient {6H{2t)). Results for collision 
angles a — 60° (o), 90° (□), and 120° ([>) are presented. The open and filled symbols represent data obtained for r = 1.0 
(fcsT varied) and ksT — 1.0 (r varied) respectively. The asterisks (*) are the theoretical predictions in the limit A/a — » 0. The 
inset shows the numerically obtained probability p2 as a function of mean free path for the different parameters considered. 
Parameters: L/a = 64, M = 5. 
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FIG. 15: Schematic diagram of contributions from 5H{2t). Solid lines, all of which have a length of two time steps, show the 
time span during which particles stay in the same collision cell and are therefore correlated. Dashed lines show the time span 
in which particles are uncorrelated, so that the molecular chaos assumption is valid. Each solid line will contribute a term 
5H{2t) and dashed lines will contribute a factor Hc{nT). 5H{nT) for n = 2, 3,4, 5 show the first four contributions. Note that 
two consecutive solid lines, i.e. as shown in 5H{At), means that particle i is correlated twice for two consecutive time steps, 
but not with the same particle. 
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FIG. 16: Schematic diagram showing ways in which two particles can be in the same shifted cell at consecutive time steps. 
The boxes with solid and dashed borders represent the shifted grids at f = r and t = 2r, respectively. S2 is the shift at t = 2r. 
The coordinate system uses the shifted frame at f = r as a reference. Two particles can be in a) the same shifted cell = 
at both t = T and t = 2r, b) cells = at t = r and = —a at t = 2t, or c) cells = at t = t and = a at t = 2t. For 
simplicity, only shifts in a;-direction are shown. 



